﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using SmartMathLibrary.Blas.Engine;

namespace SmartMathLibrary.Blas
{
    /// <summary>
    /// This class provides an interface to the Blas level 2 routines.
    /// </summary>
    [Serializable]
    public static class ManagedBlasLevel2
    {
        /// <summary>
        /// Sum of the product of a general band matrix and vector and a scaled vector. 
        /// </summary>
        /// <param name="transA">Specifies the operation to be performed on the matrix a.</param>
        /// <param name="m">The number of rows in the matrix a.</param>
        /// <param name="n">The number of columns in the matrix a.</param>
        /// <param name="kl">The lower bandwidth of the matrix a.</param>
        /// <param name="ku">The upper bandwidth of the matrix a.</param>
        /// <param name="alpha">The scalar used to multiply the matrix-vector product.</param>
        /// <param name="a">Reference to the first element in a one-dimensional array containing the elements of the matrix.</param>
        /// <param name="lda">Leading dimension of a.</param>
        /// <param name="x">A reference to a one-dimensional array containing the elements of the vector x.</param>
        /// <param name="incx">The increment for the array x.</param>
        /// <param name="beta">The scalar used to multiply y.</param>
        /// <param name="y">A reference to a one-dimensional array containing the elements of the vector y. The elements of y are overwritten with the result.</param>
        /// <param name="incy">The increment for the array y.</param>
        public static unsafe void Dgbmv(TransposeOperation transA, int m, int n, int kl, int ku, double alpha,
                                        double[] a, int lda, double[] x, int incx, double beta, double[] y, int incy)
        {
            fixed (double* pA = a)
            {
                fixed (double* pX = x)
                {
                    fixed (double* pY = y)
                    {
                        sbyte* trans;

                        if ((transA == TransposeOperation.None) || (transA == TransposeOperation.Conjugate))
                        {
                            sbyte tempuri = (sbyte) 'N';
                            trans = &tempuri;
                        }
                        else
                        {
                            sbyte tempuri = (sbyte) 'T';
                            trans = &tempuri;
                        }

                        ManagedBlasProvider.dgbmv_(trans, &m, &n, &kl, &ku, &alpha, pA, &lda, pX, &incx, &beta, pY,
                                                   &incy);
                    }
                }
            }
        }

        /// <summary>
        /// Sum of the product of a general matrix and vector and a scaled vector.
        /// </summary>
        /// <param name="transA">Specifies the operation to be performed on the matrix a.</param>
        /// <param name="m">The number of rows in the matrix a.</param>
        /// <param name="n">The number of columns in the matrix a.</param>
        /// <param name="alpha">The scalar used to multiply the matrix-vector product.</param>
        /// <param name="a">Reference to the first element in a one-dimensional array containing the elements of 
        /// the matrix.</param>
        /// <param name="lda">Leading dimension of a.</param>
        /// <param name="x">A reference to a one-dimensional array containing the elements of the vector x.</param>
        /// <param name="incx">The increment for the array x.</param>
        /// <param name="beta">The scalar used to multiply y.</param>
        /// <param name="y">A reference to a one-dimensional array containing the elements of the vector y. 
        /// The elements of y are overwritten with the result.</param>
        /// <param name="incy">The increment for the array y.</param>
        public static unsafe void Dgemv(TransposeOperation transA, int m, int n, double alpha, double[] a, int lda,
                                        double[] x, int incx, double beta, double[] y, int incy)
        {
            fixed (double* pA = a)
            {
                fixed (double* pX = x)
                {
                    fixed (double* pY = y)
                    {
                        sbyte* trans;

                        if ((transA == TransposeOperation.None) || (transA == TransposeOperation.Conjugate))
                        {
                            sbyte tempuri = (sbyte) 'N';
                            trans = &tempuri;
                        }
                        else
                        {
                            sbyte tempuri = (sbyte) 'T';
                            trans = &tempuri;
                        }

                        ManagedBlasProvider.dgemv_(trans, &m, &n, &alpha, pA, &lda, pX, &incx, &beta, pY, &incy);
                    }
                }
            }
        }

        /// <summary>
        /// Performs a rank one update of a matrix. 
        /// </summary>
        /// <param name="m">The number of rows in the matrix a.</param>
        /// <param name="n">The number of columns in the matrix a.</param>
        /// <param name="alpha">The scalar used to multiply the outer product.</param>
        /// <param name="x">A reference to a one-dimensional array containing the elements of the vector x.</param>
        /// <param name="incx">The increment for the array x.</param>
        /// <param name="y">A reference to a one-dimensional array containing the elements of the vector y. </param>
        /// <param name="incy">The increment for the array y.</param>
        /// <param name="a">Reference to the first element in a one-dimensional array containing the elements of 
        /// the matrix.</param>
        /// <param name="lda">Leading dimension of a.</param>
        public static unsafe void Dger(int m, int n, double alpha, double[] x, int incx, double[] y, int incy,
                                       double[] a, int lda)
        {
            fixed (double* pA = a)
            {
                fixed (double* pX = x)
                {
                    fixed (double* pY = y)
                    {
                        ManagedBlasProvider.dger_(&m, &n, &alpha, pX, &incx, pY, &incy, pA, &lda);
                    }
                }
            }
        }

        /// <summary>
        /// Product of a symmetric band matrix and a vector. 
        /// </summary>
        /// <param name="triangleMode">Specifies whether the matrix is an upper or lower triangular matrix.</param>
        /// <param name="n">The number of rows and columns in the matrix a.</param>
        /// <param name="k">The bandwidth of the matrix a.</param>
        /// <param name="alpha">The scalar used to multiply the matrix-vector product.</param>
        /// <param name="a">Reference to the first element in a one-dimensional array containing the elements of 
        /// the matrix.</param>
        /// <param name="lda">Leading dimension of a.</param>
        /// <param name="x">A reference to a one-dimensional array containing the elements of the vector x.</param>
        /// <param name="incx">The increment for the array x.</param>
        /// <param name="beta">The scalar used to multiply y.</param>
        /// <param name="y">A reference to a one-dimensional array containing the elements of the vector y.</param>
        /// <param name="incy">The increment for the array y. The elements of y are overwritten with the result.</param>
        public static unsafe void Dsbmv(MatrixTriangleMode triangleMode, int n, int k, double alpha, double[] a, int lda,
                                        double[] x, int incx, double beta, double[] y, int incy)
        {
            fixed (double* pA = a)
            {
                fixed (double* pX = x)
                {
                    fixed (double* pY = y)
                    {
                        sbyte* mode;

                        if (triangleMode == MatrixTriangleMode.Lower)
                        {
                            sbyte tempuri = (sbyte) 'L';
                            mode = &tempuri;
                        }
                        else
                        {
                            sbyte tempuri = (sbyte) 'U';
                            mode = &tempuri;
                        }

                        ManagedBlasProvider.dsbmv_(mode, &n, &k, &alpha, pA, &lda, pX, &incx, &beta, pY, &incy);
                    }
                }
            }
        }

        /// <summary>
        /// Performs a rank one update of a symmetric. 
        /// </summary>
        /// <param name="triangleMode">Specifies whether the matrix is an upper or lower triangular matrix.</param>
        /// <param name="n">The number of rows and columns in the matrix a.</param>
        /// <param name="alpha">The scalar used to multiply the outer product.</param>
        /// <param name="x">A reference to a one-dimensional array containing the elements of the vector x.</param>
        /// <param name="incx">The increment for the array x.</param>
        /// <param name="a">Reference to the first element in a one-dimensional array containing the elements 
        /// of the matrix.</param>
        /// <param name="lda">Leading dimension of a.</param>
        public static unsafe void Dsyr(MatrixTriangleMode triangleMode, int n, double alpha, double[] x, int incx,
                                       double[] a, int lda)
        {
            fixed (double* pA = a)
            {
                fixed (double* pX = x)
                {
                    sbyte* mode;

                    if (triangleMode == MatrixTriangleMode.Lower)
                    {
                        sbyte tempuri = (sbyte) 'L';
                        mode = &tempuri;
                    }
                    else
                    {
                        sbyte tempuri = (sbyte) 'U';
                        mode = &tempuri;
                    }

                    ManagedBlasProvider.dsyr_(mode, &n, &alpha, pX, &incx, pA, &lda);
                }
            }
        }

        /// <summary>
        /// Performs a symmetric rank two update of a symmetric matrix. 
        /// </summary>
        /// <param name="triangleMode">Specifies whether the matrix is an upper or lower triangular matrix.</param>
        /// <param name="n">The number of rows and columns in the matrix a.</param>
        /// <param name="alpha">The scalar used to multiply the outer product.</param>
        /// <param name="x">A reference to a one-dimensional array containing the elements of the vector x.</param>
        /// <param name="incx">The increment for the array x.</param>
        /// <param name="y">A reference to a one-dimensional array containing the elements of the vector y.</param>
        /// <param name="incy">The increment for the array y.</param>
        /// <param name="a">Reference to the first element in a one-dimensional array containing the elements
        /// of the matrix.</param>
        /// <param name="lda">Leading dimension of a.</param>
        public static unsafe void Dsyr2(MatrixTriangleMode triangleMode, int n, double alpha, double[] x,
                                        int incx, double[] y, int incy, double[] a, int lda)
        {
            fixed (double* pA = a)
            {
                fixed (double* pX = x)
                {
                    fixed (double* pY = y)
                    {
                        sbyte* mode;

                        if (triangleMode == MatrixTriangleMode.Lower)
                        {
                            sbyte tempuri = (sbyte) 'L';
                            mode = &tempuri;
                        }
                        else
                        {
                            sbyte tempuri = (sbyte) 'U';
                            mode = &tempuri;
                        }

                        ManagedBlasProvider.dsyr2_(mode, &n, &alpha, pX, &incx, pY, &incy, pA, &lda);
                    }
                }
            }
        }

        /// <summary>
        /// Product of a triangular band matrix and a vector. 
        /// </summary>
        /// <param name="triangleMode">Specifies whether the matrix is an upper or lower triangular matrix.</param>
        /// <param name="transA">Specifies the operation to be performed on the matrix a.</param>
        /// <param name="diag">Specifies whether or not a is unit triangular.</param>
        /// <param name="n">The number of rows and columns in the matrix a.</param>
        /// <param name="k">The bandwidth of the matrix a.</param>
        /// <param name="a">Reference to the first element in a one-dimensional array containing the elements of 
        /// the matrix.</param>
        /// <param name="lda">Leading dimension of a.</param>
        /// <param name="x">A reference to a one-dimensional array containing the elements of the vector x. The 
        /// elements of x are overwritten with the result.</param>
        /// <param name="incx">The increment for the array x.</param>
        public static unsafe void Dtbmv(MatrixTriangleMode triangleMode, TransposeOperation transA,
                                        MatrixDiagonalMode diag, int n, int k, double[] a, int lda, double[] x, int incx)
        {
            fixed (double* pA = a)
            {
                fixed (double* pX = x)
                {
                    sbyte* mode;

                    if (triangleMode == MatrixTriangleMode.Lower)
                    {
                        sbyte tempuri = (sbyte) 'L';
                        mode = &tempuri;
                    }
                    else
                    {
                        sbyte tempuri = (sbyte) 'U';
                        mode = &tempuri;
                    }

                    sbyte* trans;

                    if ((transA == TransposeOperation.None) || (transA == TransposeOperation.Conjugate))
                    {
                        sbyte tempuri = (sbyte) 'N';
                        trans = &tempuri;
                    }
                    else
                    {
                        sbyte tempuri = (sbyte) 'T';
                        trans = &tempuri;
                    }

                    sbyte* pDiag;

                    if (diag == MatrixDiagonalMode.NonUnitDiagonal)
                    {
                        sbyte tempuri = (sbyte) 'N';
                        pDiag = &tempuri;
                    }
                    else
                    {
                        sbyte tempuri = (sbyte) 'U';
                        pDiag = &tempuri;
                    }

                    ManagedBlasProvider.dtbmv_(mode, trans, pDiag, &n, &k, pA, &lda, pX, &incx);
                }
            }
        }

        /// <summary>
        /// Solves a triangular band system of equations. 
        /// </summary>
        /// <param name="triangleMode">Specifies whether the matrix is an upper or lower triangular matrix.</param>
        /// <param name="transA">Specifies the operation to be performed on the matrix a.</param>
        /// <param name="diag">Specifies whether or not a is unit triangular.</param>
        /// <param name="n">The number of rows and columns in the matrix a.</param>
        /// <param name="k">The bandwidth of the matrix a.</param>
        /// <param name="a">Reference to the first element in a one-dimensional array containing the elements of
        /// the matrix.</param>
        /// <param name="lda">Leading dimension of a.</param>
        /// <param name="x">The x.</param>
        /// <param name="incx">The increment for the array x.</param>
        public static unsafe void Dtbsv(MatrixTriangleMode triangleMode, TransposeOperation transA,
                                        MatrixDiagonalMode diag, int n, int k, double[] a, int lda, double[] x, int incx)
        {
            fixed (double* pA = a)
            {
                fixed (double* pX = x)
                {
                    sbyte* mode;

                    if (triangleMode == MatrixTriangleMode.Lower)
                    {
                        sbyte tempuri = (sbyte) 'L';
                        mode = &tempuri;
                    }
                    else
                    {
                        sbyte tempuri = (sbyte) 'U';
                        mode = &tempuri;
                    }

                    sbyte* trans;

                    if ((transA == TransposeOperation.None) || (transA == TransposeOperation.Conjugate))
                    {
                        sbyte tempuri = (sbyte) 'N';
                        trans = &tempuri;
                    }
                    else
                    {
                        sbyte tempuri = (sbyte) 'T';
                        trans = &tempuri;
                    }

                    sbyte* pDiag;

                    if (diag == MatrixDiagonalMode.NonUnitDiagonal)
                    {
                        sbyte tempuri = (sbyte) 'N';
                        pDiag = &tempuri;
                    }
                    else
                    {
                        sbyte tempuri = (sbyte) 'U';
                        pDiag = &tempuri;
                    }

                    ManagedBlasProvider.dtbsv_(mode, trans, pDiag, &n, &k, pA, &lda, pX, &incx);
                }
            }
        }

        /// <summary>
        /// Product of a triangular matrix and a vector. 
        /// </summary>
        /// <param name="triangleMode">Specifies whether the matrix is an upper or lower triangular matrix.</param>
        /// <param name="transA">Specifies the operation to be performed on the matrix a.</param>
        /// <param name="diag">Specifies whether or not a is unit triangular.</param>
        /// <param name="n">The number of rows and columns in the matrix a.</param>
        /// <param name="a">Reference to the first element in a one-dimensional array containing the elements of
        /// the matrix.</param>
        /// <param name="lda">Leading dimension of a.</param>
        /// <param name="x">The x.</param>
        /// <param name="incx">The increment for the array x.</param>
        public static unsafe void Dtrmv(MatrixTriangleMode triangleMode, TransposeOperation transA,
                                        MatrixDiagonalMode diag,
                                        int n, double[] a, int lda, double[] x, int incx)
        {
            fixed (double* pA = a)
            {
                fixed (double* pX = x)
                {
                    sbyte* mode;

                    if (triangleMode == MatrixTriangleMode.Lower)
                    {
                        sbyte tempuri = (sbyte) 'L';
                        mode = &tempuri;
                    }
                    else
                    {
                        sbyte tempuri = (sbyte) 'U';
                        mode = &tempuri;
                    }

                    sbyte* trans;

                    if ((transA == TransposeOperation.None) || (transA == TransposeOperation.Conjugate))
                    {
                        sbyte tempuri = (sbyte) 'N';
                        trans = &tempuri;
                    }
                    else
                    {
                        sbyte tempuri = (sbyte) 'T';
                        trans = &tempuri;
                    }

                    sbyte* pDiag;

                    if (diag == MatrixDiagonalMode.NonUnitDiagonal)
                    {
                        sbyte tempuri = (sbyte) 'N';
                        pDiag = &tempuri;
                    }
                    else
                    {
                        sbyte tempuri = (sbyte) 'U';
                        pDiag = &tempuri;
                    }

                    ManagedBlasProvider.dtrmv_(mode, trans, pDiag, &n, pA, &lda, pX, &incx);
                }
            }
        }

        /// <summary>
        /// Solves a triangular system of equations. 
        /// </summary>
        /// <param name="triangleMode">Specifies whether the matrix is an upper or lower triangular matrix.</param>
        /// <param name="transA">Specifies the operation to be performed on the matrix a.</param>
        /// <param name="diag">Specifies whether or not a is unit triangular.</param>
        /// <param name="n">The number of rows and columns in the matrix a.</param>
        /// <param name="a">Reference to the first element in a one-dimensional array containing the elements of
        /// the matrix.</param>
        /// <param name="lda">Leading dimension of a.</param>
        /// <param name="x">The x.</param>
        /// <param name="incx">The increment for the array x.</param>
        public static unsafe void Dtrsv(MatrixTriangleMode triangleMode, TransposeOperation transA,
                                        MatrixDiagonalMode diag,
                                        int n, double[] a, int lda, double[] x, int incx)
        {
            fixed (double* pA = a)
            {
                fixed (double* pX = x)
                {
                    sbyte* mode;

                    if (triangleMode == MatrixTriangleMode.Lower)
                    {
                        sbyte tempuri = (sbyte) 'L';
                        mode = &tempuri;
                    }
                    else
                    {
                        sbyte tempuri = (sbyte) 'U';
                        mode = &tempuri;
                    }

                    sbyte* trans;

                    if ((transA == TransposeOperation.None) || (transA == TransposeOperation.Conjugate))
                    {
                        sbyte tempuri = (sbyte) 'N';
                        trans = &tempuri;
                    }
                    else
                    {
                        sbyte tempuri = (sbyte) 'T';
                        trans = &tempuri;
                    }

                    sbyte* pDiag;

                    if (diag == MatrixDiagonalMode.NonUnitDiagonal)
                    {
                        sbyte tempuri = (sbyte) 'N';
                        pDiag = &tempuri;
                    }
                    else
                    {
                        sbyte tempuri = (sbyte) 'U';
                        pDiag = &tempuri;
                    }

                    ManagedBlasProvider.dtrsv_(mode, trans, pDiag, &n, pA, &lda, pX, &incx);
                }
            }
        }
    }
}